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EIGENVECTOR-BASED CENTRALITY MEASURES FOR 
TEMPORAL NETWORKS * * * § 
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Abstract. Numerous centrality measures have been developed to quantify the importances of 
nodes in time-independent networks, and many of them can be expressed as the leading eigenvector 
of some matrix. With the increasing availability of network data that changes in time, it is important 
to extend such eigenvector-based centrality measures to time-dependent networks. In this paper, we 
introduce a principled generalization of network centrality measures that is valid for any eigenvector- 
based centrality. We consider a temporal network with N nodes as a sequence of T layers that 
describe the network during different time windows, and we couple centrality matrices for the layers 
into a supra-centrality matrix of size NT x NT whose dominant eigenvector gives the centrality of 
each node i at each time t. We refer to this eigenvector and its components as a joint centrality , as 
it reflects the importances of both the node i and the time layer t. We also introduce the concepts 
of marginal and conditional centralities, which facilitate the study of centrality trajectories over 
time. We find that the strength of coupling between layers is important for determining multiscale 
properties of centrality, such as localization phenomena and the time scale of centrality changes. 
In the strong-coupling regime, we derive expressions for time-averaged centralities , which are given 
by the zeroth-order terms of a singular perturbation expansion. We also study first-order terms to 
obtain first-order-mover scores, which concisely describe the magnitude of nodes’ centrality changes 
over time. As examples, we apply our method to three empirical temporal networks: the United 
States Ph.D. exchange in mathematics, costarring relationships among top-billed actors during the 
Golden Age of Hollywood, and citations of decisions from the United States Supreme Court. 
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1. Introduction. The study of centrality 94 —that is, determining the impor¬ 
tances of different nodes, edges and other structures in a network —has widespread 
applications in the identification and ranking of important agents (or interactions) 
in a network. These applications include ranking sports teams or individual ath¬ 
letes 14 16 110 , the identification of influential people 170 , critical infrastructures 
that are susceptible to congestion 53 59 , impactful United States Supreme Court 


cases j4lj|42 ,79 , genetic and protein targets 65 , impactful scientific journals |7 , and 
much more. An especially important type of centrality measure are ones that arise as 
a solution of an eigenvalue problem 46 , with nodes’ importances given by the entries 


of the dominant eigenvector of a so-called centrality matrix. Prominent examples in¬ 
clude eigenvector centrality 9], PageRank 12 95] (w hich provides the mathematical 
foundation to the Web search engine Google 46,78]), hub and authority (i.e., HITS) 
73 , and Eigenfactor |7 . However, despite the fact that real-world networks 


scores 


change with time 60-62 


most methods for centrality (and node rankings that are 
derived from them) have been restricted to time-independent networks. Extending 
such ideas to time-dependent networks (i.e., so-called temporal networks) has recently 
become a very active research area [T][ 33 l[ 5 T][ 52 l[ 7 T][Ml[ 90 l[ 9 ^ 1 fl 07 l| 108 Ull 6 yil 7 || 127 || 128 l . 

In the present paper, we develop a generalization of eigenvector-based centrali¬ 
ties to ordered multilayer networks such as temporal networks. Akin to multilayer 
modularity 92 , a critical underpinning of our approach is that we treat inter-layer 
and intra-layer connections as fundamentally different 26 72 . For the purpose of 


generalizing centrality, we implement this idea by coupling centrality matrices for 
temporal layers in a supra-centrality matrix. For a temporal network with N nodes 
and T time layers, we obtain a supra-centrality matrix of size NT x NT. To use this 
matrix, we represent the temporal network as a sequence of network layers, which is 
appropriate for discrete-time temporal networks and continuous-time networks whose 
edges have been binned to produce a discrete-time temporal network. Importantly, 
our methodology is independent of the particular choice of centrality matrix, so it 
can be used, for example, with ordinary eigenvector centrality 9], hub and authority 


scores 


73 , PageRank 95 , or any other centralities that are given by components of 


the dominant eigenvector of a matrix. 

The dominant eigenvector of a supra-centrality matrix characterizes the joint 
centrality of each node-layer pair [i, t )—that is, the centrality of node i at time step 
t —and thus reflects the importances of both node i and layer t. We also introduce 
the concepts of marginal centrality and conditional centrality , which allow one to (1) 
study the decoupled centrality of just the nodes (or just the time layers) and (2) 
study a node’s centrality with respect to other nodes’ centralities at a particular time 
t (i.e., the centralities are conditional on a particular time layer). These notions make 
it possible to develop a broad description for studying nodes’ centrality trajectories 
across time. Although we develop this formalism for temporal networks, we note 
that our approach is also applicable to multiplex and general multilayer networks, 
which are two additional scenarios in which the generalization of centrality measures 
is important. (See, e.g., 51 112 113 for multiplex centralities and 130,114 for general 
multilayer centralities.) 

Similar to the construction of supra-adjacency matrices 


49 72 92 (see Sec. 2.2 


for definition), we couple nearest-neighbor temporal layers using inter-layer edges of 
weight w, which leads to a family of centrality measures that are parameterized by 
ui. The parameter oj controls the coupling of each node’s eigenvector-based centrality 
through time and can be used to tune the extent to which a node’s centrality changes 
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Fig. 1.1. As one case study for our generalization of eigenvector-based centralities to temporal 
networks, we examine a multilayer temporal network that we construct from data obtained from 
the Mathematics Genealogy Project (MGP) he- (See the Supplementary Material or nn to 
download the temporal network data.) We study a temporal network with T = 65 time layers 
(corresponding to the years 19j6-2010), in which a given edge j —»■ i in time layer t signifies the 
number of graduating doctoral students in year t at university i who later supervise a graduating 
doctoral student at university j. The nodes’ sizes and colors indicate what we call “time-averaged 
centrality” (see Sec. 0- which we calculate for the type of centrality matrix known as an authority 
matrix 1 7 if . We do not show self-edges. Because visualizing a temporal network is difficult \25\ , we 
show the network that we obtain by aggregating the adjacency matrices across time layers. Studying 
such an aggregated network neglects the time-ordered structure of a temporal network. In our paper, 
we study trajectories of node centralities (which represent importances) over time. See Sec. \ 5.1 for 
additional discussion of the MGP network. (We created this image using the software Gephi 1 3'.) 


over time. The limiting cases ui —> 0 + and uj —> oo are particularly interesting, as 
the former represents the regime of complete decoupling of the layers and the latter 
represents a regime of dominating coupling of the layers. As part of the present paper, 
we conduct a perturbative analysis for the u> —> oo limit. (For notational convenience, 
we study e —> 0 + for e = 1/w.) This allows us to derive principled expressions for 
the nodes’ time-averaged, centralities (given by the zeroth-order expansion) and first- 
order-mover scores , which are derived from the first-order expansion. Time-averaged 
centrality ranks nodes so that their centralities are constant across time, and first- 
order-mover scores rank nodes according to the extent to which their centralities 
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change in time. The computation of both time-averaged centralities and first-order- 
mover scores can be very efficient, because they only require the numerical solution 
of linear-algebraic problems of size N x N, which is ordinarily much smaller than the 
full supra-centrality matrix of size NT x NT. Moreover, our perturbative approach 
also alleviates the need to demand a particular choice of intra-layer coupling weight co. 
Note that our methodology is a retrospective analysis and complements alternative 
temporal generalizations that take a causal approach. 

The remainder of our paper is organized as follows. In Sec. [2j we provide further 
background information. In Sec. [3] we present our generalization of eigenvector-based 
centralities for temporal networks. In Sec. [4] we derive principled expressions for time- 
averaged centralities and first-order-mover scores based on a singular perturbation 
expansion. In Sec. [5j we illustrate our approach using three examples from empirical 
data: Ph.D. exchange using data from the Mathematics Genealogy Project [101 


(see Fig. top billing in the Golden Age of Hollywood using the Internet Movie 
Database (IMDb) [63], and citations of United States Supreme Court decisions 40 
64 . We conclude in Sec. [6] and provide further details of our perturbation expansion 
in the appendix. 


2. Background Information and a Naive Approach. In Sec. |2.1[ we pro¬ 
vide additional background information. In Sec. |2.2[ we discuss a naive approach 
for temporal eigenvector centrality that motivates our approach, and introduces our 
mathematical notation and terminology. 

2.1. Background Information. The analysis of large networks is ubiquitous in 
science, engineering, medicine, and numerous other areas 94 . In the social sciences, 


for example, the abundance of data that describes the social behavior of individuals 
in academia 13 15 19 93 97 101 , show business 111 , politics 8 41 42 79 100] , 


and just about every other arena offers exciting avenues for the quantitative study of 
social systems. For these and many other applications, it is important to develop (and 
improve) mathematical techniques to extract concise and intuitive information from 
large network data. From the interdisciplinary pursuit of what is now often called 
network science, we know—from theory, computation, and data analysis—that many 
network properties (e.g., degree heterogeneity, local clustering, community structure, 
and others 94 ) have significant effects on dynamical processes on networks (e.g., 
information dissemination and disease spreading) [20] [99] [T 25 ] . Although the vast 


majority of research in network science has focused on time-independent networks, 
increased effort in recent years has aimed to generalize network analyses to “temporal 
networks” |4|6^[62][69], in which network entities and/or interactions change in time. 
In the study of time-independent networks, numerous centrality measures have 


been developed to try to quantify the relative importances of nodes 

94 

125 . There 

are seemingly as many centrality measures as applications [6| 9 

10, 

34 l 

43 

67 , and 


different types of centrality are appropriate for different situations. Importantly, one 
can construct centrality measures not only from the direct consideration of network 
structure but also based on studying an appropriate dynamical process on a net¬ 


work 45 95 105 . In our paper, we focus on eigenvector-based centralities. Although 


eigenvectors can obviously be used in different ways to introduce different notions of 
centrality, we reserve the term eigenvector-based centrality to refer only to centrality 
measures in which the nodes’ centralities are given by the entries of the dominant 
eigenvector of a matrix, which we call a centrality matrix. Centrality matrices in¬ 
clude a network’s adjacency matrix A (which indicates which nodes share a common 
edge) and various functions of A 6|[35], such as the hub and authority matrices 73 
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and PageRank matrices 46 78 95 . Eigenvector-based centralities reflect a network’s 


global structure and are often preferable to other types of centralities, because there 
are a large number of computationally efficient algorithms for computing spectral 
properties of matrices (e.g., the power method for computing the dominant eigenvec¬ 
tor 48 ). An eigenvector-based centrality is also a key component to the function of 


major Web search engines, including Google 46 73 95 


Although the study of centrality measures in time-independent networks has been 
insightful for numerous applications, most networks are time-dependent, and it is 
important to generalize centrality measures for temporal networks. This has been a 
very active area of research in recent years. Many approaches consider a temporal 
network as a sequence of network layers, and initial studies examined centralities 
of such uncoupled individual network layers 11 . It was found, for example, that 


in the absence of coupling between network layers, centrality scores can fluctuate 
significantly from one time step to the next due to stochasticity in the appearance 
and disappearance of edges. Rankings can of course change over time (either slowly or 
rapidly). To give an example that is directly relevant to one of our case studies, it was 
demonstrated recently that uncoupled university rankings can fluctuate significantly 
from year to year 115 , which we consider to be a suboptimal description for the 


temporal evolution of departmental prestige. 

The most popular temporal extensions of centrality measures involve consider¬ 
ing notions of centrality based on paths in a time-independent network and gener¬ 
alizing them using time-respecting paths in a temporal network 74 75 . We point 


out temporal extensions to PageRank that have been introduced to counteract age 
bias 85 124 128 , capture temporal changes for the external interests of nodes via 

and study random walks taking place on temporal net- 


108 


dynamic teleportation 

works [108|. We also highlight temporally extended versions of betweenness cen¬ 
trality TTT171 Trf 127 


closeness centrality 71 96 


win/lose centrality 90], communicability (33|50]52 , Katz centrality 51 , and coverage 
centrality 116 . These research efforts have extended centrality measures in differ- 


117 127 , Bonacich centrality 80 


ent ways, which one can relate partly to the fact that one can define time-respecting 
paths in different ways. A time-respecting path can allow one, several, or an unlim¬ 
ited number of edge traversals during a particular time step. Moreover, the length 
of a time-respecting path—which one can use to provide a notion of “distance” be¬ 
tween nodes—can also be measured in different ways 127 . For instance, the length 


of a time-respecting path can describe the number of edges that are traversed by the 
path, latency between the initialization and termination times of the path, or a com¬ 
bination of such ideas. For temporal networks, even the notion of a “shortest path” 
lacks an unambiguous definition, which subsequently can also introduce ambiguity 
into dynamical processes and network measures derived from them. 

One common feature of existing generalizations of centralities for temporal net¬ 
works is that they illustrate the importance of studying an entire temporal network, 
as opposed to aggregating temporal layers into a single (time-independent) network 
or analyzing the time layers in isolation from one another 60|61|. Specifically, study¬ 
ing a layer-aggregated network prevents one from studying centrality trajectories (i.e., 
how centrality changes over time), and studying the time layers in isolation does not 
account for the temporal orderings of edges, which can be crucial for determining 
centralities in a temporal setting 053,381150 51 84 . To provide additional context 


we 


highlight that dynamical processes can behave vastly differently on temporal versus 
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layer-aggregated networks j]] For example, a random walk—a process on which many 
eigenvector-based centralities rely—on a temporal network is affected fundamentally 
by the temporal ordering and time scale of the appearances and disappearances of 
edges 57 58 60-62 66 . Rankings, such as eigenvector-based centralities, that are 
derived from such dynamics are, in turn, affected fundamentally by the temporal 
structure of the networks, and aggregation (as well as isolation) can lead to mislead¬ 
ing or even simply wrong results. Additionally, if one starts with a Markovian process 
on a temporal network and then aggregates the network, then in general one does not 
obtain a Markovian process 62 , so fundamental (and often desirable) properties of a 


dynamical process can be destroyed as a byproduct of neglecting a network’s inherent 
temporal structure. 

Finally, although extending centrality to temporal networks has become a very 
active research area, one should expect different generalizations of eigenvector-based 
centralities to be appropriate for different applications. Existing papers have not 
always been clear about the modeling assumptions and tradeoffs of their approaches. 

2.2. Naive Generalization of Eigenvector Centrality for Temporal Net¬ 
works. Before we present our primary approach in Sec. [3], it is instructive to consider 
one possible way to generalize eigenvector-based centralities. This example motivates 
our approach and introduces mathematical notation and terminology. Importantly, 
it is naive in that it does not treat intra-layer edges and inter-layer edges as distinct 
types of edges, which causes problems when the network layers are strongly coupled. 

We use a multilayer representation of networks and seek to identify the most 
central nodes of a temporal network with N distinct nodes (i.e., vertices or actors) 
across T time layers. We specify the network edges with a node-by-node-by-time 
(N x N x T) adjacency tensor in which nonzero elements indicate the presence 
and weight of the edge from node i to node i in time layer t. That is, the adjacency 
matrix at time t is given by AW See Table 


2.1 


for a summary of our mathematical 
notation. We refer to node i in layer t as a node-layer pair (i, t ) and node i (regardless 
of layer) as a physical node. We are interested particularly in understanding the 
physical nodes’ centrality trajectories through time. 

It is tempting to reshape a network’s associated adjacency tensor into an NT x NT 


1 As discussed in, e.g., |26|7 2 and several references therein—and more recently in |27| —a similar 
issue arises more generally in multilayer networks, and one must also take into account the effects of 
inter-layer edges (which are fundamentally different from intra-layer edges) when defining a dynamical 
process on a multilayer network. 


Table 2.1 

A summary of our mathematical notation. 


Typeface 

Class 

Dimensionality 

M 

matrix 

NT x NT 

M 

matrix 

NxN 

M 

matrix 

TxT 

V 

vector 

NT x 1 

V 

vector 

N x 1 

V 

vector 

T x 1 

Mij 

scalar 

1 

Vi 

scalar 

1 
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supra-adjacency matrix 


- A d) 

col 

0 

col 

At 2 ) 

col 

0 

col 

A (3) 


( 2 . 1 ) 


which represents a collection of both the temporal network edges (i.e., intra-layer 
edges) and the “identity edges” (which are inter-layer edges) that couple the node¬ 
layer pairs {(z,f)} for the same physical node i across the T network layers. The 
identity edges of weight co attempt to weight the persistence of a physical node through 
time by enforcing an identification with itself at consecutive times [5]. When there 
are inter-layer edges only between different instances of the same physical node, a 
multilayer network is said to exhibit diagonal coupling, and the use of a constant 
co across all such inter-layer edges is sometimes known as layer coupling 72 . We 


restrict our attention to nonnegative inter-layer coupling to > 0. (One can consider 
w < 0 to drive negative coupling between layers, but we do not examine such values 
in our applications.) One can construe w as a parameter to tune interactions between 
network layers [2,5 91 . In the limit co —> 0 + , the layers become uncoupled; in the 


limit w —>■ oo, the layers are so strongly coupled that inter-layer weights dominate the 
intra-layer connections. 

We also restrict ourselves to nearest-neighbor coupling of temporal layers, as we 
place the identity inter-layer edges only between node-layer pairs, (z, t ) and ( i , t ± 1), 
that are adjacent in time (where the t = 0 and t = T layers have inter-layer edges 
to one other layer rather than two). This results in the block structure in Eq. (2.1). 
Equivalently, we write 


A = diag f A W.A^l + Ad hain ) (8) col , 


( 2 . 2 ) 


where 0 denotes the Kronecker product and A.) chain ) is the T xT adjacency matrix of 
an undirected chain , or “bucket brigade,” network whose T nodes are each adjacent to 
their nearest neighbors along an undirected chain. In this bucket brigade, = 1 

for j = i ± 1 and A^ halr ^ = 0 otherwise. Although one can choose inter-layer coupling 
matrices other than A) cham ) for the inter-layer couplings 24 72 (and much of our 


approach can be generalized to other choices of coupling), we restrict our attention 
to nearest-neighbor coupling of layers. 

It is also tempting to directly apply a standard eigenvector-based centrality cal¬ 
culation to the supra-adjacency matrix A by treating it just like any other adjacency 
matrix despite its special structure. However, such an approach neglects to respect 
the fundamental distinction between intra-layer edges and inter-layer edges that arises 
from the block-diagonal structure of A. That is, in such an approach, one treats 
the inter-layer couplings (i.e., identity arcs) just like any other edge. In general, 
however, one needs to be careful when studying a temporal network using a supra- 
adjacency matrix formalism because many basic network properties—some of which 
carry strong implications about a time-independent network (e.g., its spectrum, con¬ 
nectedness properties, and so on)—do not carry over naturally without modification 
to the supra-adjacency matrix. This issue was discussed for multilayer networks more 


generally in Refs. 21 26 72 and more recently in Ref. 27 
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As a concrete example, let’s examine hub and authority centralities (i.e., HITS 


731) for a directed temporal network using the supra-adjacency matrix in Eq. (2.11 by 


simply inserting it in place of a time-independent adjacency matrix in the standard 
formulas. In other words, we define the hub and authority matrices as AA T and A T A, 
respectively. At a glance, by noting that the inter-layer couplings are undirected but 
that the intra-layer edges are directed, we already see that it is not clear whether 
standard interpretations of hub and authority rankings are still sensible. Nevertheless, 
one can try this approach for computing generalized hub and authority scores as the 
dominant eigenvectors of the symmetric matrices AA T and A T A. The simplicity of 
this approach makes it pleasing (and tempting), and the two symmetric matrices do 
have a block structure. However, in contrast to A (whose blocks on and off of the 
main diagonal encode intra-layer and inter-layer edges, respectively), the blocks in the 
matrices AA T and A T A no longer separate neatly into describing only a single type 
of edge (i.e., inter-layer versus intra-layer edges). 

The problem with this construction becomes particularly clear in the limit of 
strong inter-layer coupling (i.e., as w -> oo), for which A « w(A/ chain ) g I). Because 
A (chain) ig symme tric, it follows that AA T « A T A « w 2( A (chain)( A (chai„))T ^ j). 

Unfortunately, it is useless to compute hub and authority scores of an undirected 
chain. Specifically, the corresponding hub/authority centrality matrix (whose domi¬ 
nant eigenvector gives the hub/authority scores) of the undirected chain becomes 


(chain) (chain) 'yT _ (chain) ^T^(chain) _ 


1010 
0 2 0 1 


1 0 
0 1 


2 0 
0 2 


0 

1 0 


0 


(2.3) 

revealing that the hub/authority scores of the even-indexed and odd-indexed nodes 
decouple from each other. The resulting matrix is no longer irreducible, which can 
lead to nonuniqueness of dominant eigenvectors and/or can also cause the entries of 
a dominant eigenvector to be identically 0 for a large number of nodesj^j Both issues 
are detrimental if one wants to rank nodes based on some notion of importance. For 
example, for large values of w, we observe oscillations and numerical instabilities when 
attempting to generalize hub and authority centralities in this way. 


2 By inspection, the matrix in Eq. I 2 T 3 J is not irreducible, so we cannot apply the Perron- 
Frobenius theorem for irreducible nonnegative matrices [87] . There exist many variations of Perron- 
Frobenius theory, including ones that are applicable to reducible matrices (e.g., see Sec. 8.3.1 of [87] ), 
which one can use to study phenomena that arise in the absence of irreducibility. In our case, we find 
two types of scenarios, depending on whether N is odd or even. For even N, the largest eigenvalue of 
A(chain) ^(cham)^T h as a multiplicity of two and a corresponding two-dimensional eigenspace that 
is spanned by vectors in which either the even-indexed or odd-indexed entries are 0. Consequently, 
there is not a unique dominant eigenvector, so there is not a unique ranking of nodes. For odd N , 
there is a single dominant eigenvalue; however, its eigenvector has entries that are identically 0 for 
even-indexed nodes, so only half of the nodes are ranked in a nontrivial way. 
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3. Temporal Coupling of Eigenvector-Based Centralities. In this section, 
we present a mathematical formalism for eigenvector-based centralities in temporal 
networks that treats inter-layer and intra-layer edges as distinct types of edges and 
ensures appropriate behavior for all to > 0. Similar to prior investigations using mul¬ 
tilayer representations of temporal networks 72 92 , we seek to develop an approach 


that involves neither a heuristic averaging of centralities from individual layers nor 
invokes the centrality for a single network obtained from the aggregation of network 
layers (e.g., summing the network edges across time). 

The remainder of this section is organized as follows. In Sec. |3.1[ we present 
our methodology for temporal eigenvector-based centrality in terms of the dominant 
eigenvector of a supra-centrality matrix. In Sec. |3.2| we introduce the concepts of joint, 
marginal, and conditional centrality, and we use them to study decoupled centralities 
of nodes and layers based on the centralities of node-layer pairs. In Sec. |3.3[ we 
illustrate these concepts for an example synthetic network. 


3.1. Inter-Layer Coupling of Centrality Matrices. To avoid the problems 
that arise from ignoring the distinction between inter-layer edges and intra-layer edges, 
we define a somewhat more nuanced generalization of eigenvector-based centralities. 
To preserve the special role of inter-layer edges, we directly couple the matrices that 
define the eigenvector-based centrality measure within each temporal layer (e.g., or¬ 
dinary adjacency matrices for eigenvector centrality). That is, one can cast any 
eigenvector-based centrality in terms of some matrix C that is a function of the 
adjacency matrix A. For example, hub and authority scores are the leading eigenvec¬ 
tors of the matrices AA r and A 1 A, respectively (using the convention that elements 
Ajj indicate i —> j edges). Letting C- 1 denote the centrality matrix for layer t, we 
couple these centrality matrices with inter-layer couplings of strength w in a (rescaled) 
supra-centrality matrix 


C(e) 


' eCW I 0 

I eC( 2 ) I 

0 I eC^ 


(3.1) 


We are defining the supra-centrality matrix using a scaling factor e = 1/ui, because 
it gives convenient mathematical notation for our forthcoming singular perturbation 
analysis. However, because it is more intuitive (and standard) to describe layers as 
being coupled together by weight w, we refer to the limit e —> 0 + (i.e., oj —> oo) as the 
strong-coupling regime. 

One can interpret the parameter e > 0 as a tuning parameter that controls how 
strongly a given physical node’s centrality is coupled to itself between neighboring time 
layers. (See the related discussions in j2j : 5 in the context of multilayer community 
structure.) That is, the intuition for a specified eigenvector-based centrality proceeds 
within each individual layer as in the associated centrality’s original definition, and 
the additional inter-layer coupling introduces contributions to centrality from the 
network structure in neighboring layers. Of particular interest are the limits in which 
e->-oo (i.e., decoupling of layers) and e —> 0 + (i.e., a particular notion of order¬ 
preserving aggregation). See the related discussions in [102 103] . We expect the 
e —> 0 + limit to yield principled time-averaged centralities of nodes. Note that such a 
notion reflects the layers having an intrinsic temporal ordering and should in general 
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yield different results from calculating the centralities of summed adjacency layers (i.e., 
directly summing the corresponding entries in these matrices) or from an unweighted 
averaging of centralities across otherwise uncoupled layers. 

The notation of Eq. (3.1) assumes implicitly that every physical node i appears in 
every time layer t (where i £ {1,..., N} and is{l,...,T}). Although this notation 
is consistent with that used in the development of multilayer modularity 


92 


it is 

important to call attention to practical issues in treating situations in which some 
physical nodes do not appear across all layers 21 26 72 . When defining multilayer 


modularity, there are no difficulties with removing nodes from the layers in which they 
do not appear as long as the correct identity inter-layer edges are coded appropriately. 
(Indeed, see the U.S. Senate roll-call voting example of 92].) In contrast, for reasons 
that will become clearer as we develop our singular perturbation analysis in the strong¬ 
coupling limit (see Sec. [4]), temporal generalization of eigenvector-based centralities 
in this limit requires that all physical nodes are taken into account across all layers, 
even when they do not appear in a layer. In other words, one must account for each 
physical node i as a “ghost node” in any layer in which i does not appear in the data. 
The ghost node is adjacent to its counterparts in neighboring layers via inter-layer 
edges, and it therefore maintains connectivity to the full multilayer network; however, 
it does not have any intra-layer edges, as such edges are “forbidden.” 

Before continuing, we briefly comment on assumptions that we make about the 
supra-centrality matrix C(e). In the construction of eigenvector-based centralities 
for time-independent networks, it is typically assumed that a given centrality matrix 
is nonnegative and irreducible 6, 35}[46jj73,78, 95 . Similarly, we assume that C(e) 
is nonnegative and irreducible for any e > 0. Our motivation is that the Perron- 
Frobenius theorem for nonnegative matrices [87 ensures that the largest (positive) 
eigenvalue has multiplicity one and that its corresponding eigenvector is nonnegative 
and unique, which are both beneficial properties for ranking nodes based on a notion 
of importance. Similar to the case of time-independent networks, one can guarantee 
that the matrix C(e) is both nonnegative and irreducible by placing simple constraints 
on the properties of the temporal network. For example, consider the matrix C(e) 
and its associated network -that is, the network in which every nonzero entry of C(e) 
corresponds to an edge from some node i to some other node j, and the weight is 
given by the value of the entry. (In practice, one can use such an approach to study 
any matrix if one interprets its nonzero entries in terms of a network.) It follows that 
C(e) is irreducible and nonnegative if this associated network is strongly connected. A 
sufficient (but not necessary) condition to assure this is that all centrality matrices 
are nonnegative and the aggregated matrix CW itself has an associated network 
that is strongly connected. For example, when computing eigenvector centrality for an 
undirected temporal network (i.e., C® = A^), this constraint implies that > 
0 and that the aggregation J2t °f the adjacency matrices yields an adjacency 
matrix with an associated strongly-connected network]^] In general, however, the 


3 There exist both stronger and weaker versions of such a relation between network structure 
and the dominant eigenspace of the matrices that are associated with a network. If a network is 
strongly connected, then the largest (positive) eigenvalue of the matrix has a multiplicity of one, 
and its corresponding eigenvector is guaranteed to be unique and strictly positive. If a network is 
weakly connected and if all nodes are contained in the union of the largest in-, out-, and strongly- 
connected components, then the matrix has a largest (positive) eigenvalue with a multiplicity of one, 
and its corresponding eigenvector is both unique and nonnegative. In other cases, the eigenvector 
corresponding to the largest eigenvalue may or may not be unique. If there are edges with negative 
values, the eigenvector corresponding to the largest eigenvalue may not be nonnegative. 
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irreducibility of C(e) depends on that of the centrality matrices (i.e., Y2t C (t ^), rather 
than on the adjacency matrices. 

Our assumption that C(e) be irreducible places restrictions both on the set {C^} 
and on our choice for how to couple the layers. We choose to couple each time 
layer t to both layers t + 1 and t — 1 (when present), so our temporal extension of 
eigenvector centrality is non-causal; that is, the centralities are coupled both forward 
and backward in time. In principle, one could choose other strategies for coupling the 
layers. Note, however, that a causal strategy in the absence of other features (e.g., 
one could examine different types of “teleportation” [46]) forces C(e) to be reducible, 
because the centralities of past network layers cannot depend on future layers. 

3.2. Joint, Marginal and Conditional Centrality for Multilayer Net¬ 
works. We study the dominant eigenvector v(e) of C(e), with corresponding (and 
largest positive) eigenvalue A max (e) [i.e, C(e)v(e) = A max (e)v(e)]. The entries of the 
dominant eigenvector give the centralities of each node-layer pair (z,f); this repre¬ 
sents the centrality of physical node i at time t. The dominant eigenvector v(e) of a 
supra-centrality matrix in Eq. ( |3.1| ) gives the centrality of node-layer pairs. That is, 
the eigenvalue entry ujv(t-i)+i(c) indicates the centrality of node i at time t. Such 
a joint centrality, whether given by w(e) or any other centrality for node-layer pairs, 
reflects information about the importances of both the nodes and the layers. We 
develop a simple formalism to decouple these centralities. For concreteness, we use 
v(e), but this approach can be applied to any centrality measure of node-layer pairs 
in a multilayer (e.g., temporal) network including those not based on eigenvectors. 

Our approach is inspired by multivariate statistics: we define “joint”, “marginal”, 
and “conditional” centralities. Joint centrality describes the importances of node-layer 
pairs, marginal centrality describes the uncoupled centrality of either nodes or layers, 
and conditional centrality describes the importance of a node-layer pair as compared 
to, for example, other node-layer pairs in that same layer. 

To proceed, it is convenient to map the vector w(e), which has length NT, to an 
N x T matrix W, which we define entry-wise by 


Wa — v N (t-i)+i(e) ■ (3-2) 

The scalar W it gives the joint centrality of the node-layer pair (i, t); that is, it indicates 
the centrality of node i at time t. We define the marginal node centrality (MNC) Xi 
and marginal layer centrality (MLC) yt by 

x t = - yt = Yl w * • ( 3 - 3 ) 

t i 

The values {xt} and {yt} indicate the importances of nodes and layers, respectively, 
for a particular choice of e. Although we use the summation to compute marginal 
node and layer centralities, one can also consider other aggregation methods. We 
define the conditional centrality of node-layer pair (*,t), conditioned on layer t, by 

Zit = Wit/yt ■ (3.4) 

The scalar Z lt indicates the importance of physical node i relative to other physical 
nodes in layer t. For some applications, it can be beneficial to similarly study the 
conditional centrality of layers conditional on a given node, but we do not explore 
this notion in the present paper. For a given node i £ {1, ...,1V} and time t G 
{1,..., T}, the sets {Wit} and {Zu} of centrality values indicate trajectories for how 
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(a) Temporal Network (b) 


Centralities for e = 0.5 
Time Layer 





Fig. 3.1. Eigenvector centralities for an example undirected temporal network, with N = 4 
nodes and T = 3 layers, in which we use the layers’ adjacency matrices as the centrality matrices. 

(a) A temporal network with intra-layer edges (black lines) and inter-layer identity edges (gray lines). 

(b) For e = 0.5, we indicate the joint node-layer centralities { IT,, } (white boxes), which correspond 
to the entries in v(Y) and represent the centralities of each physical node in each layer (there are 
N X T such centralities). We refer to the marginal centralities {x, } and {yt} (shaded boxes) that 
one obtains by summing rows and columns, respectively, as the “marginal node centralities" (MNC) 
and “marginal layer centralities’’ (MLC). (c) MNC and MLC versus the coupling parameter e. For 
small e, the dominant time layer is the central layer (2’); for large e, the dominant layer is layer 
1 ’, which contains the centrality matrix with the dominant (i.e., largest positive) eigenvalue. The 
vertical dashed line indicates the value e = 0.5 that we use in panels (b) and (d). (d) We examine 
the nodes’ centrality trajectories by plotting the joint node-layer centralities (upper subpanel) and 
conditional node-layer centralities (lower subpanel) versus time. 


the importance of physical node i changes through time. We interpret conditional 
centrality trajectories as follows: for a given physical node i, we study a sequence of 
centralities in which the fth term indicates a centrality that is relative to centralities 
of node-layer pairs at time t. This contrasts with the joint node-layer centralities; 
because ||wr(e) || 2 = 1, joint node-layer centralities indicate a notion of importance 
that is relative to all node-layers pairs. 

3.3. Synthetic Temporal Network Example. We illustrate the concepts of 
joint, marginal, and conditional centrality with a toy example. The network, which 
we show in Fig. |3.1[a), consists of N = 4 physical nodes and T = 3 time layers. 
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We compute eigenvector centralities of this (undirected) temporal network by setting 
each layer’s centrality matrix to be its adjacency matrix and solving the dominant- 
eigenvalue equation, C(e)v(e) = A max (e)v(e), for several choices of e. In Fig. |3.1[ b), 
we summarize the centrality measures for e = 0.5 in a matrix. The entry in row i and 
column t gives Wn. The MNCs {a;,} and the MLCs {y t } are in the shaded boxes, and 
they indicate the relative centralities of the physical nodes and layers, respectively, 
for the chosen value of e. 

In Fig. |3.1[ c), we plot the dependence of the MNC (upper subpanel) and MLC 
(lower subpanel) on the coupling strength e. For small e, the dominant time layer 
is layer 2’; for large e, the dominant layer is layer 1’. This is unsurprising: when 
considering the centrality matrices of the layers in isolation, the dominant eigenvalue 
of the centrality matrix of layer 1’ is larger than that for the other layers. The choice 
of e is very important, and centralities can depend discontinuously on e because of 
eigenvalue crossings. In this example, we obtain three qualitatively different regimes: 
(i) a strong-coupling regime in which the centralities are similar to what we expect in 
the limit e —» 0 + ; (ii) a weak-coupling regime in which the centralities behave similarly 
to what they do in the e —> oo limit; and (iii) an intermediate-coupling regime in which 
the centralities transition between these two limiting cases. (Compare this result to 
the phase-transition phenomena for graph Laplacians of multilayer networks discussed 
in Refs. [102,103] .) We also explore these regimes for our case study with the MGP 
network in Sec. EU 

In Fig. injd), we plot joint node-layer centralities (upper subpanel) and the con¬ 
ditional node-layer centralities (lower subpanel) that correspond to the four physical 
nodes. We show results for e = 0.5. Note that the conditional node-layer centralities 
for physical node 4 increase over time, whereas they decrease for physical node 1. 
This is sensible for the temporal network in Fig. 3.1 a): at time t = 1’, physical node 


1 has the largest degree and physical node 4 has the smallest degree; however, at time 
t = 3’, physical node 1 has the smallest degree and physical node 4 has the largest 
degree. 

4. Singular Perturbation in the Strong-Coupling Limit. Given the joint 
node-layer centralities and conditional node-layer centralities that correspond to a 
physical node i , it is possible to define a notion of “time-averaged centrality” by 
summing one of these centralities across the layers (e.g., the prior yields the MNC). 
However, it is not clear which is preferable, and these centralities are sensitive to the 
value of e. Alternatively, one can define a time-averaged centrality by studying the 
limit e —> 0 + . In this limit, the conditional node-layer centrality of every physical 
node i is constant across the time layers. 

Examining centralities as e —> 0 + provides a principled approach for calculat¬ 
ing time-averaged centralities. However, the supra-centrality matrix C(e) given by 
Eq. (3.1) becomes singular at e = 0, which complicates the consideration of this limit. 


The intra-layer connectivity is completely eliminated, and the network decomposes 
into N connected components. That is, the matrix is no longer irreducible, so the 
Perron-Frobenius theorem no longer holds. Indeed, at e = 0, the dominant eigenvalue 
has an TV-dimensional eigenspace. In contrast, for e > 0, the dominant eigenvalue has 
a single eigenvector. 

To overcome this issue, we derive a singular perturbation expansion in the limit 
e —► 0 + . In Sec. 4.1 we further explore the singularity that arises in the strong¬ 


coupling limit. In Secs. |4~2| and [4~3| we give zerotli-order and first-order perturbation 
expansions, which lead to principled expressions for time-averaged centralities and 
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first-order-mover scores, respectively. We give higher-order expansions in an appendix. 
In Sec. |4.4[ we summarize our procedure and discuss the computational complexity 
of computing time-averaged centralities and first-order-mover scores. 


4.1. Singularity at Infinite Inter-Layer Coupling. In this section, we de¬ 
velop a perturbation analysis of the dominant eigenspace (i.e., the eigenspace of the 
largest positive eigenvalue) of C(e) [see Eq. (3.1)] in the limit e —► 0 + . To allow this 
analysis to have more broad applicability, we do a perturbation expansion using the 
following general form of coupling of block matrices: 


1(e) = B + eG , 


(4.1) 


where G = diagfM^ 1 ),..., M*^ T 1], B = A <g> I, and the T x T matrix A (recall Table 
2 .1) encodes the inter-layer coupling, where entry A tt / indicates how layer t is coupled 


to layer t' . We recover the supra-centrality matrix C(e) in Eq. (3.1) by using nearest- 
neighbor layer coupling, A = ^4( chain ) ) and centrality matrices along the diagonal 
(i.e., Mw = C^). Additionally, similar to our assumptions for Eq. (3.1), we assume 
that M(e) is nonnegative and irreducible for any e > 0. These assumptions hold as 
long as the summation of matrices (y~] f M^) and inter-layer coupling matrix A each 
correspond to a strongly connected network (see footnote 3 in Sec. |3.2[ ). 

We begin by studying the dominant eigenspace for the matrix M(e) at e = 0. We 
thus consider the matrix 


= A®I. 


(4.2) 


and we will show that the dominant eigenspace of M(0) is iV-dimensional, implying 
that one cannot obtain the unique dominant eigenvector of M(e) for the e —>■ 0 + limit 
simply by setting e = 0. Instead, one needs to do a singular perturbation analysis. 
To facilitate our discussion, we define an NT x NT stride permutation matrix P i|48 
with entries Pm = 1 for l = \k/N ] +T [(fc— 1) mod IV] and Pm = 0 otherwise. (Recall 
that [•] denotes the ceiling function.) Note that P permutes node-layer indices, so we 
can easily go back and forth between ordering the node-layer pairs first by time and 
then by physical node index, or vice versa (i.e., ordering them first by physical node 
index and then by time). The main benefit of defining P is that one can express . 
in the form M(0) = P (I <g> A) P T 28 , where 


I® A = 


A 0 
0 A 

0 0 


(4.3) 


Clearly, I ® A decouples into N identical eigenvalue equations for the inter-layer 
coupling matrix A. Because P is a unitary matrix, this decoupling also occurs for 
M(0). Specifically, the NT eigenvalues of I (g> A are given by the T eigenvalues of 
A , where each eigenvalue has multiplicity N and a corresponding iV-dimensional 
eigenspace spanned by the vectors that can be constructed using the eigenvectors of 
A (i.e., with appended 0 values in appropriate coordinates). 

We explain this construction in more detail for the dominant eigenspace. Let 
v denote the largest positive eigenvalue of A, and let it = [ui,... ,mt] T denote its 
corresponding eigenvector. Because of the block-diagonal structure of Eq. (4.3), the 
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largest eigenvalue Ao of 1(g) A is Ao = v, and its corresponding eigenspace is spanned by 
the N eigenvectors {uj, where o, = [0 T ,..., 0 T , it T , 0 r ,... , 0 T ] 5 . That is, the ith 
block of iij is given by u , and all of the other blocks are vectors of zeros. Consequently, 
one can obtain the N dominant eigenvectors of M(0) = A ® I using the permutations 
{Pui}. That is, they have the general form ^WayPuj, where the constants {ai} 
must satisfy ]Tb af = 1 to ensure that the vector is normalized. Because there does 
not exist a unique dominant eigenvector of M(0) since its dominant eigenspace is N 
dimensional, we need to develop a singular perturbation analysis to obtain a unique 
solution for the dominant eigenvector of M(e) in the e —> 0 + limit. 

When the network layers are coupled by an undirected chain network, the inter¬ 
layer coupling matrix is given by A = ^ 4 ( cham ) ^ which has N eigenvalues and eigen¬ 


vectors given by 82 


^(chain) = 2 cog 


(chain)_ 


T+ 1 


y/ln 


sin 


T17T 

T+ 1 


, sin 


2mr 
T + 1 


, sm 


Tmr 

T+ 1 


i T 


(4.4) 

(4.5) 


where n £ {1,..., N} and the normalization constant is y n = Y^t= i sin 2 [mrt/(T + 1)]. 
Setting n= 1 gives the dominant eigenvalue and its corresponding eigenvector.. 

4.2. Zeroth-Order Expansion and Time-Averaged Centrality. In this 
section, we study the zeroth-order expansion of the dominant eigenvector w(e) in 
the limit e —» 0 + . As we now show, the conditional node-layer centralities {Zn} 
corresponding to a given physical node i become constant across time in this limit. 


(Recall our definitions of joint, marginal and conditional centralities in Sec. 3.2 ) For 
each physical node, we refer to this limiting value as its time-averaged centrality. By 
examining the first-order expansion of v(e) in the limit e —> 0 + , we show in Eq. (4.121 


that one can obtain the time-averaged centralities as the eigenvector components cor¬ 
responding to the largest eigenvalue of a matrix of size N x N. 

We consider the dominant-eigenvalue equation 


Amax(eMe) = M(e)v(e) = Bw(e) + eGv(e). 


(4.6) 


We expand A max (e) and v(e) for small e by writing A 
v(e) = wq + ewi + ■ • • to obtain /cth-order approximations: A 


(e) — Ao T eAi T * * * and 

e) ~ Ej=o eiA i and 

v ( e ) ~ Y^ k j=o e ^j■ W e use superscripts to indicate powers of e in the terms in the 
expansion, and we use subscripts for the terms that are multiplied by the power of e. 
Note that Ao and w a respectively indicate the dominant eigenvalue and corresponding 
eigenvector in the limit e —> 0 + . Successive terms in these expansions represent higher- 
order derivatives of A max (e) and w(e), and each term assumes appropriate smoothness 
of these functions. 


Our strategy is to develop consistent solutions to Eq. (4.6 1 for increasing values 
of k. Starting with the first-order approximation, we substitute A max (e) ss Ao + eAi 
and w(e) ss v 0 + evi into Eq. (4.6) and collect the zeroth-order and first-order terms 
in e to obtain 


(A 0 I —B)v 0 = 0, (4.7) 

(AqI — B) = ((& — Ail) "vq , ( d -$) 
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where I is the NT x NT identity matrix. Equation (|4.7|) is exactly the system that we 


studied in Sec. 4.1 [see Eq. (4.1 1 with e = 0], where we found that the operator Aol — B 
is singular and has an iV-dimensional null space. (This is the dominant eigenspace of 


.) We also found that Eq. (4.71 has a general solution of the form 


An — 


= £ 


ay-Puj , 


(4.9) 


where {ai} are constants that satisfy the constraint that wy has magnitude 1 (i.e., 
J2i a i = !)• We defined Uj just below Eq. (4.3). 


To find the set {a,} of unique constants that determine w 0 , we seek a solvability 
condition in the first-order terms. Using the fact that the null space of AqI — B is 


span(Pui,... ,Pujv) for any physical node i, it follows that 
and left-multiplying Eq. (4.8) by (Pui) T leads to 


ufP T (GW/o = Aiui 


TmT, 


huiP (A 0 I - B) vi = 0, 


(4.10) 


Using the solution of wq in Eq. (4.9), we obtain 


ajTufP T GPuj = Ai ajTufP T PTUj = Ai on , 


(4.11) 


because P T P = PP T = I and ujiij = 5ij, where Sij is the Kronecker delta. Letting 
a = [ai,..., aiv] T , Eq. (4.11 1 corresponds to an iV-dimensional eigenvalue equation, 


xWq = Aia, 

where the matrix X* 1 ! has elements 


(t) 2 
’u t . 


(4.12) 

(4.13) 


Our assumption that M(e) is nonnegative and irreducible for any e > 0 ensures that 
X« is also nonne gat ive and irreducible. By the Perron-Frobenius theorem for non¬ 


negative matrices 87 , the largest positive eigenvalue Ai of X* 1 ^ has a multiplicity of 


one, and its eigenvector a is unique and has nonnegative entries. (See Sec. 3.1 and 
footnote[3j) We normalize the solution a to Eq. (4.12) by = 1 an< ^ substitute 

the normalized solution into Eq. (4.9) to obtain the zeroth-order term v q. 


When the layers are coupled by an undirected chain and the block matrices are 
the layers’ centrality matrices (i.e., = C^), we obtain 




(t) ■ 2 

ij sm 


7T t 

T + 1 


(4.14) 


where 71 = iSin 2 ( 7 rt/(T + 1 )) is the normalization constant for the dominant 
eigenvector tt( chain ) given by n = 1 in Eq. (4.51. In this case, recall that the vector 
wo is the dominant eigenvector of C(e) in the limit e —> 0 + and gives the joint node¬ 
layer centralities in this limit. By inspection, we see that the elements of wq are 
a i sin( 7 rt/(T + 1)) for node-layer pair ( i,t ). Because these correspond to the limiting 
e —>■ 0 + entries of vector "v(e), they are independent of e. At the same time, conditional 
centrality of node-layer pair (*,f) is ai (up to a normalization constant), independent 
of the layer t. That is, the conditional node centrality trajectories become constant 
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across time in the limit e —► 0 + . These values of {cti} arise naturally from our 
perturbative expansion in the supra-centrality framework, independently of the value 
of e. By contrast, recall that the marginal node centralities (MNCs) reflect averaging 
the joint centralities across time layers for a specific choice of e. Accordingly, we 
hereafter refer to the entry at in the vector a as the time-averaged, node centrality 
of physical node i. Because our approach can also be applied to multilayer networks 
that are not necessarily temporal, we call the layer-averaged node centrality for 
these situations. 


4.3. First-Order Expansion and First-Order-Mover Scores. In this sec¬ 


tion, we show that the first-order expansion of Eq. (4.61 leads to a linear system [see 
Eq. (4.24)], which we solve to obtain a measurement of the variation over time of each 
physical node’s centrality trajectory [see Eq. (4.26)]. Specifically, as one increases 
e above 0 + , the first-order expansion, (which accounts for first derivatives with re¬ 
spect to e), captures the dominant changes in centrality trajectories for small values 
of e. (In the appendix, we derive expressions for higher-order terms that account for 
higher-order derivatives with respect to e.) 

In Sec. |4.2[ we derived closed-form expressions for Ao and and an eigenvalue 
equation satisfied by Ai. We now solve for vj to complete our first-order approxi¬ 
mation. For notational convenience, we define Lo = Aol — B and Li = G — Ail, so 
Eq. (4.8) becomes Lq^i = Li^ 0 - Letting Lj denote the Moore-Penrose pseudoinverse 


of Lo, we write 


ii = LqLiwo + /SjPuij = LqGvq + /3j-Pu 


(4.15) 


We simplify the first term in Eq. (4.15) using Li = G — Ail and Vo = Ylj o^Puij an d 
noting that each vector Pu 7 lies in the null space of each of the matrices Lo and Lj. 
The second term in Eq. (4.15) accounts for the projection of vi onto the null space 
of Lo, where the constants /?, = (Pui) T vi indicate the projections onto the spanning 
vectors of the null space. To ensure numerical stability and computational efficiency 
in practice, we calculate Lj using the identity 


28 


Lj = (A 0 I - A ® I) f = (A 0 I — A) + G I. 


(4.16) 


Note that Lq depends only on the inter-layer coupling matrix A (e.g., for nearest- 
neighbor-in-time coupling, A = A.( chain )), which one can compute and save in memory 
prior to analyzing network data. 

Just as we examined first-order terms to solve for constants {a,;}, we now seek 


a solvability condition in the second-order terms to determine {/3i} in Eq. (4.15). 
Substituting A max (e) = Ao + eAi + e 2 A 2 and v(e) = vo + evi + e 2 ^2 into Eq. (4.6) and 
collecting the second-order terms yields 


L 0 tv 2 = Liwi - A2-VQ ■ 


(4.17) 


Similar to before, we left-multiply Eq. (4.17) by (Pu;) T and require both sides to be 
identically 0 to obtain 


A 2 ufP 7 Vo 
Using a* = urf P T ^o and /3, 


: UlfP T Ll'Vl = lufP T G'Vl — AlinfP T ’Vl. 

(4.18) 

: TufP 7 ^!, it then follows that 


A 2 g^ + Ai /3i = ufP T G’vi. 

(4.19) 
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Substituting the expressions for wo from Eq. (4.9) and wi from Eq. (4.151 into 
Eq. (4.19) then yields 

\ 20 ti + Ai/3j = ^ ajufP T GLQ®PiUj + ^ /3,-uifP T GPiu :) -. (4.20) 

0 j 

After some rearranging, we obtain 

(X« - ArI)/3= (A 2 I-X( 2 ))a, (4.21) 

where the matrix X^ 1 ) was defined by Eq. (4.13), and the elements of the matrix 'SS 2 ' 1 
are 


a1 2) = u7 1 p t gl);gpiu 1 . 


(4.22) 


Recalling that we determined a as the solution of X^cc = Aicc such that JU a 2 = 1, 
we left-multiply Eq. (4.21) by a T to obtain 


A 2 =a T X (2 )c*. 


(4.23) 


We thereby obtain 


(3 = (X^ - Ai^^Aal - X^)a + ba , (4.24) 


where the constant b = a T (3 describes the (possibly nonzero) projection of (3 onto 
the null space of (X^ — Ail) [see Eq. (|4.12[ )]. 

We now show that b = 0 in Eq. ( |4.24[ ), by virtue of the requirement that the 
eigenvector obtained at first-order has a norm of 1. That is, we require that 1 = 
||v 0 + evi|| 2 = ||v 0 || 2 + 2e(v 0 ,vi) + 0(e 2 ). However, ||tv 0 || 2 = 1, so (v 0 ,vi) = 0, 
where we use the notation (•, •) to denote the dot product between inputs. Using the 
definitions of "Vq and ¥ 1 , we see that 


0 = (vo, "Vi) = ( ^2 ajPuj, LqG'Vo + /%IPtu-j 

\ i o 

= 51 a 3^3 = b > 


(4.25) 


because the vectors {Piu ; y} are orthonormal and lie in the null space of Lo- (In other 
words, (Pu,;) T Pui :) - = 5ij and (Puj) T LQ = 0.) 

In practice, we solve Eq. ( |4.24[ ) using a linear solver (see, e.g., [48]) rather than 
the pseudoinverse to avoid computing the inverse of (X^ — Ail). We then ensure that 
the solution is orthogonal to a by projecting it onto the subspace that is orthogonal 
to a. 


One can substitute the solution f3 to Eq. (4.24) with 5 = 0 into Eq. (4.15) to 


obtain the first-order term wi in the expansion for v(e). This first-order term, which 
yields the strongest temporal variation of the conditional centralities at small e, is a 
concise representation of temporal changes in centrality. There are multiple possible 
ways to use vi to quantify the role of physical node i across the T layers. We define a 
measure rrii that equals the square root of the sum of the squares of the entries in vi 
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that correspond to physical node i. Specifically, we define the first-order-mover score 
mi > 0 of physical node i by 

to? = vfPI i P T ^i 

= P? + (lJg^o) T PI j P T U <g ^o 

T 2 

= Pi + ^ ([LjGv 0 ],: + t(jv_i)^ , (4.26) 


where [■]* denotes the itli entry in a vector and I; = diag[0,..., 0,J,0,..., 0] is a 
matrix of size NT x NT that contains all 0 entries except for the zth block, which 
is an identity matrix I of size T x T. In other words, we measure the variation 
of w-| with respect to a physical node i by examining the 2-norm of the entries in 
vj that correspond to the node-layer pairs (i,t) that are relevant to physical node 
i (i.e., entries j such that j = N(t — 1) + i, with t = 1,... ,T). In principle, one 
can also use a different vector norm or a heuristic method for aggregating centrality. 
Our choice has the virtue that it is mathematically consistent with our definition for 


the time-averaged centralities {at} [see Eq. (4.121]. Specifically, af = v(fPIjP T vo. 
Therefore, one can naturally extend our approach for quantifying the contribution 


of the first-order correction given by Eq. (4.26) to higher-order corrections. We also 


note that first-order-mover scores rank the nodes according to the magnitudes of their 
corresponding entries in vi- Therefore, the associated centralities w(e) can either 
increase or decrease over time. One can easily check whether there is an increase or 
a decrease by examining the corresponding entries of v(e). 

4.4. Procedure for Computing Time-Averaged Centrality and First- 
Order-Mover Scores. We summarize our procedure for computing time-averaged 
centrality and first-order-mover scores: 

1. Construct the matrix X^ 1 - 1 using Eq. (4.14): 


^= 7l - 1 V.d?sin 2 


t 


nt 

T + 1 


When layers are coupled by a layer-adjacency matrix A (which is not nec¬ 
essarily an undirected chain 24.( cham )), it follows that = uJ’P^’GPuj, 

where G = diag[C^\ ..., C (2 )] and P and tu* are defined just before and 
after, respectively, Eq. (14.31). 


2. Solve for the time-averaged centralities {c^} using Eq. (4.12): 

X (1) a = Aia. 

3. Construct the matrix X( 2 ) using Eq. (4.22): 


Xi 2) = u^GliGPiu 


3 » 


where Lq = (A 0 I — A)^ ® I. 

4. Solve for (3 in Eq. ( |4.21 1: 

(X (1) - AiI)/3= (A 2 I-X (2) )a. 


where A 2 = a T X( 2 ^o:. 
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5. Solve for the first-order-mover scores {mt} using Eq. (4.261: 


2 

m i = Pi + ([^o < G'v 0 ] i+t ( JV _i)J , 

t =i 

where "Vo = JT a^-Piiij. 

We comment briefly on the computational costs of the above procedure. The 


supra-centrality matrix [see Eq. (3.1)], whose dominant eigenvector gives the joint 
node-layer centralities, has size NT x NT , and that can be problematic for large net¬ 
works with many time layers (i.e., when T 1). The time-averaged node centralities 


are given by the solution to Eq. (4.12), which is a dominant eigenvalue problem for a 


matrix of size N x TV. To examine which physical nodes have centralities that change 
significantly over time, we examine the first-order-mover scores given by Eq. (|4.26 ); 
this requires one to solve the iV-dimensional linear system given by Eq. (4.24). Be¬ 


cause Lq, (G, and vq are known prior to solving Eq. (4.24), one can directly compute 


the second term in Eq. (4.26). For sparse networks [i.e., those in which the number 


of edges at a given time is 0(iV)], the matrices that we have discussed in this section 
are typically also sparse. One can thus solve Eqs. (4.12), (4.14), (4.21), (4.22), and 


(4.26) efficiently using data structures that are designed for sparse matrices, includ- 
iterative methods 


ing direct methods [23 
network structures (e.g., 


109 


and methods designed for particular 
nested dissection for planar networks [8l]). In particular, 


the power method for computing a dominant eigenvalue and eigenvector of a sparse 
matrix reduces the per-iteration complexity from 0(N 2 ) to O(M), where M is the 
number of nonzero entries in the sparse matrix. However, the actual scaling can be 
much larger, because the number of iterations required for convergence depends on 
the gap between the largest and second-largest eigenvalues. 

5. Case Studies with Empirical Network Data. In this section, we examine 
temporal centrality in case studies with three sets of empirical data: the Mathematics 


Genealogy Project (MGP; see Sec. 5.1) of Ph.D. receipt in the mathematical sciences 
in U.S. universities, top billing in the Golden Age of Hollywood (GAH; see Sec. 5.2), 
and citations of U.S. Supreme Court decisions (SCD; see Sec. 5.3). We have posted 
Matlab software at 118 that implements our calculations and can be used to re¬ 


produce the results of this section. 

Most of our calculations for these examples use a temporal generalization of hub 
and authority scores [73], which are particularly appropriate for directed networks 
(such as our three examples). We also note that hub and authority scores have 
been used previously to examine time-independent faculty-hiring networks 39,93 and 
Supreme Court decisions |39][4l]. To illustrate a comparison with another choice of 
centrality, we also study a temporal generalization of PageRank for the MGP network. 

We compute hub and authority scores independently using two supra-centrality 
matrices, A^[A^] 7 and [A^J^AW For a single-layer network, it is possible to 
simultaneously compute hub and authority scores by studying the single system 
0 . In general, however, a supra-centrality matrix that uses this alternative 
formulation yields different time-dependent centralities from ones based on indepen¬ 
dent computations of hub and authority scores. For each centrality that we compute, 
we also examine the induced ranking of nodes in which the highest-ranked nodes (i.e., 
those with ranks 1, 2, and so on) correspond to the largest centralities, whereas the 
lowest-ranked nodes correspond to the smallest centralities. 
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5.1. Doctoral Degree Exchange in the Mathematics Genealogy Project 
(MGP). Our first case study uses a network that encodes the exchange of mathemati¬ 
cians (and other mathematical scientists) who have obtained a Ph.D. (or equivalent 
doctoral degree) between universities to study the academic prestige of those univer¬ 
sities. We study data provided by the Mathematics Genealogy Project (MGP) [101 


which collects information for mathematicians (and members of related fields who are 
listed in the MGP) with doctorates. For each mathematical scientist, the informa¬ 
tion includes graduation year, his/her official academic advisor(s), the degree-granting 
university, and a list of his/her students who have also obtained doctoral degrees. A 
subset of the present authors previously utilized this data to approximate the flow of 
doctorates between universities—that is, a person graduates from one university and 
is then hired at a second university—and quantified the resulting hub and authority 
scores for the total flow during a specified time period as a candidate measure of these 
universities’ relative mathematical prestige 93 . Moreover, hub and authority scores 


have been used previously to study Ph.D. exchange networks in other disciplines 39 


See 

55 

for 

see 

44 

83] 


for a comparison of the hiring market for different academic disciplines, and 
for other analyses and visualization of data from the MGP. See 76 for 


an application of PageRank centrality to ranking world universities using data from 
Wikipedia. 

It is well-documented that graduates typically obtain faculty positions at univer¬ 
sities that are either comparable to or less prestigious than the one from which they 
graduate 


13 15 19 29 39 97 


and we study university prestige as indicated by the 
exchange between universities of mathematicians with doctoral degrees. We general¬ 
ize a previous study of the MGP data in 


93 


by keeping the year that each faculty 
member graduated with his/her Ph.D. degree. We focus on the years 1946-2010, 
which includes all post-World War II information available in the data set]/] This 
yields T = 65 time layers, and we restrict our attention to a set of N = 227 U.S. 
universities that were connected during this period. To construct the network, we 
create directed intra-layer edges i —► j at time t to represent a doctoral degree in the 
MGP data awarded to a mathematical scientist from university j in year t who later 
advised at least one student at university i. Therefore, to contribute a directed edge, 
the mathematician must have at least one student in the MGP data. We weight edges 
to indicate the number of doctorates from university j in year t awarded to faculty 
who later advise students at university i. Our construction aligns edge directions to 
be opposite to that of the flow of people, so a node with large in-degree (i.e., with 
many graduates who later advise students elsewhere) is considered both an academic 
authority as well as an authority with respect to HITS centrality 731. See Fig. |l.l| for 
a visualization of this network; due to the difficulty of visualizing temporal networks 
25 , we show a network corresponding to the aggregation Y2 t of the adjacency 
matrices. Although one could define the multilayer network in more intricate ways 
(e.g., by normalizing edge weights using the number of graduates) and examine how 
the results vary for different choices, we wish to keep the present manuscript focused 
on introducing and demonstrating our temporal generalization of eigenvector-based 
centralities. Therefore, we leave such detailed analyses for future work. We make the 
MGP temporal network available as a Supplementary Material and online at 


119 


4 The data set was provided to us in 2009, although it includes information up to 2010. The year 
2006 is the last year in which a Ph.D. degree was awarded to someone who was subsequently a Ph.D. 
advisor in the data, so it is also the last year in which intra-layers edges are present. Additionally, 
we decided to be optimistic and include Ph.D. degrees that were projected for the year 2010. 
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Table 5.1 

Top time-averaged centralities and first-order-mover scores for our temporal generalization [see 
Eq. ( |3. ![ >/ of authority scores | ? ;[ for U.S. universities in the Mathematics Genealogy Project \l 0 1L 
We show our results using two different orderings: (1) according to the top time-averaged centralities 
and (2) according to the top first-mover scores. 


Top Time-Averaged Centralities Top First-Order-Mover Scores 


Rank 

University 

Cti 

mi 

University 

rm 


1 

MIT 

0.6685 

688.62 

MIT 

688.62 

0.6685 

2 

UC Berkeley 

0.2722 

299.06 

UC Berkeley 

299.07 

0.2722 

3 

Stanford 

0.2295 

241.71 

Princeton 

248.72 

0.1803 

4 

Princeton 

0.1803 

248.71 

Stanford 

241.71 

0.2295 

5 

UIUC 

0.1645 

74.30 

Georgia Tech 

189.34 

0.0960 

6 

Cornell 

0.1642 

180.50 

U Maryland 

186.65 

0.1278 

7 

Harvard 

0.1628 

185.34 

Harvard 

185.34 

0.1628 

8 

U Washington 

0.1590 

81.22 

CUNY 

182.59 

0.0466 

9 

U Michigan 

0.1521 

86.50 

Cornell 

180.50 

0.1642 

10 

UCLA 

0.1456 

152.77 

Yale 

159.11 

0.0816 


5.1.1. MGP: Centrality in the Strong-Coupling Regime. We begin by 
identifying the universities that have the largest time-averaged authority centrali¬ 
ties {ai}, which we obtain from the dominant eigenvector of the matrix [see 

Eq. (4.14)] using C^ = (A ^) 2 A^ 1 73 . For notational convenience, we use t to de¬ 
note the graduation year rather than the time layer. For example, we use A^ 1946 ) to 
denote the network adjacency matrix for time layer 1 (i.e., year 1946). We summarize 
these authority values in Table |5.1[ and we note that the most central universities 
according to this measure are all widely-accepted top-tier programs in mathematics. 
The time-averaged authorities identify the four most central mathematics universities 
for this time period as MIT, UC Berkeley, Stanford, and Princeton. Although the 
results in Table 5.1 are interesting, time-averaged centrality (by definition) does not 
provide information about temporal trajectories of the universities’ authorities, and 
this is the type of idea that we seek to explore. We thus calculate the first-order-mover 
scores {rrii} from Eq. (4.26), and we list the universities with the top first-order-mover 

Note the similarity between the two lists; that 


scores in the right column of Table 5.1 
is, universities with largest ai tend to also have large m 


In Fig. 5.1 we show further results for prestige (as revealed by Ph.D. exchange). 
In Fig. |5.1[ a), we plot university ranking according to {trii} versus its ranking accord¬ 
ing to {«i}. Note in the bottom left corner that MIT is ranked first for both quan¬ 
tities, and that in general there is a strong linear correlation between rank according 
to (Xi and rank according to mi. Intuitively, this suggests that shifts in centrality 
include a natural effect that is related directly to the centrality score itself. (In other 
words, large centrality values tend to also include large fluctuations, whereas small 
centrality values typically have only small fluctuations.) Deviations from the observed 
nearly-linear relation indicate universities whose centrality trajectory exhibits larger 
variations over time, and it is worthwhile to look at these universities in more detail 
for potentially interesting insights. For example, the universities with high rank ac¬ 
cording to rrii (i.e., large m*) but comparatively low rank according to oti (i.e., small 
cti) include Georgia Tech and CUNY, and it is known that Georgia Tech’s mathe¬ 
matics department transitioned from a primarily teaching-oriented department to a 
much more research-oriented department with a newly restructured doctoral degree 
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time-averaged centrality rank 



Centrality trajectories for e — 10 

MIT 

UC Berkeley 
Princeton 
Stanford 
Georgia Tech^ 
CUNY 


1945 


1965 1985 

time ( t ) 


2005 


Fig. 5.1. University rankings in the Mathematics Genealogy Project (MGP) HE according 
to time-averaged centralities and first-order-mover scores. We give results for our temporal gener¬ 
alization [see Eq. ( |3.1| )7 of authority scores \73\ . (a) We plot the first-order-mover ranking of nodes 
(i.e., ranked according to {mi}) versus the time-averaged ranking of nodes (i.e., ranked according 
to {cx-i}). As shown in the inset, nodes with large time-averaged rank tend to also have large first- 
order-mover rank (e.g., MIT ranks first in both). However, there are nodes that have a much higher 
first-order-mover rank than time-averaged centrality rank (e.g., Georgia Tech and CUNY). In panel 
(a), we show a magnification of the gray box in the subpanel, (b) We plot the conditional node-layer 
centralities (i.e., the centrality of node-layer pairs normalized by the centrality of each time layer) to 
study the universities ’ centrality trajectories over time. We show results for some of the top ranked 
first-order-movers. Most of these top first-order-movers are also top time-averaged authorities (e.g., 
MIT). In contrast, Georgia Tech and CUNY rank in the top six of the first-order-movers ranking, 
but they are in the lower reaches of the top 4.0 for the time-averaged ranking. As expected, this 
ranking difference reflects the fact that their centrality trajectories exhibit a significant change over 
time. Georgia Tech rises in rank when t E [1965,1985], whereas CUNY’s rank drops during this 
time period. 


program starting in the late 1970s 32 


In Fig. 5.1 b), we plot the conditional authority centralities at e = 10 4 of univer¬ 


sities versus time for six of the universities with the largest first-order-mover scores 
rrii. This includes the four universities with the top time-averaged centralities, as 
well as Georgia Tech and CUNY (which do not have highly-ranked time-averaged 
centralities). As we expect, the conditional centralities for Georgia Tech and CUNY 
change drastically over time, whereas the trajectories for the others remain relatively 
constant. 

5.1.2. MGP: Some Properties of Authority Centrality. As we showed in 


Sec. 3.3 for a synthetic network, the choice of inter-layer coupling strength e strongly 
affects the temporal behavior of a node’s centrality trajectory. We expect to observe 
three qualitative regimes: (i) the strong-coupling (e —k 0 + ) regime that we studied 
in Sec. [4j (ii) a weak-coupling (e —> oo) regime; and (iii) an intermediate-coupling 
regime, in which the centralities behave differently than expected for either the strong¬ 
coupling or weak-coupling regimes. In Fig. |5.l[ b), we show results for e = 10~ 4 , and 
we observe that the universities tend to have slowly-varying centrality trajectories. 
However, the choice of e should depend both on the application and on the question 
of interest. As we are about to illustrate, it is important to consider what values of e 
are appropriate. 

In Fig. |5.2[ w e study centrality trajectories for Georgia Tech for various choices 
for e. In Fig. |5.2[ a), we show the conditional node-layer centralities for Georgia Tech 
versus time t. Recall that the conditional node-layer centrality indicates the centrality 
of node-layer pair (i,t) with respect to all node-layer pairs (j, t) at time t. We also 
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(b) Georgia Tech joint centrality 



Fig. 5.2. Centrality trajectories for Georgia Tech illustrate that one can construe e as a tuning 
parameter that controls how much centrality can vary between neighboring time layers, (a) To study 
the trajectory of university authorities over time, we examine the conditional node-layer centralities. 
For sufficiently small e, we observe a steady increase in ranking with time for Georgia Tech. Varying 
e changes the coupling strength between temporal layers. If e is too large (e.g., e > 0.01J, then the 
coupling between layers is so weak that the conditional node-layer centrality of the two node-layer 
pairs at times t and t + 1 are no longer similar in value. As e —¥ 0“*", the conditional node-layer 
centrality approaches the stationary, time-average ranking given by oti (horizontal dashed line), but 
we still observe significant variation even for e = 0.0001. (b) We plot the joint node-layer centrality 
of the T node-layer pairs that correspond to Georgia Tech across the time layers for several values of 
e. For small e, the joint node-layer centralities are determined by the chain of identity edges (which 
leads to the sinusoidal dependence given by Eq. ( |4.5| with n = 1 and magnitude OLi). For large 
values of e, the node-layer pairs in time layer t = 1982 dominate the joint node-layer centralities, 
so all universities have their highest values in this layer. 


show the value of oti (rescaled for normalization), which gives the conditional node¬ 
layer centrality of Georgia Tech in the limit e —> 0 + . For small but nonzero e (e.g., 
e = 10 -3 ), note that we obtain a similar trajectory as for e = 0 + . For example, the 
trajectory varies slowly over time, so the conditional node-layer centrality of Georgia 
Tech at times t and t + 1 are approximately equal for all t. However, as we increase 
e, we lose the slow temporal variation over time. For example, when e > 10 1 , the 
conditional centrality of Georgia Tech at times t and i + 1 are typically very dissimilar, 
which appears to be a consistent property of conditional node-layer centralities in the 
limit e — > oo. It is our believe that the highly volatile rankings for large e do not 
appropriately describe the dynamics of department prestige 115] ; this observation has 
motivated us to focus on the small e (i.e., strong coupling) regime in this paper. The 
limiting cases e —> 0 + and e —> oo, respectively, do a good job of describing regimes 
with very small and very large e, but the intermediate (“transitional”) regime between 
these extremes is not straightforward to interpret. Even the boundaries between the 
two extreme qualitative regimes are not clear and are open to interpretation. 

In Fig. |5.2[ b), we plot the joint node-layer centralities for Georgia Tech for various 
values of e. Recall that the joint node-layer centrality of the node-layer pair (i. t) 
reflects information about both the physical node i and the time layer t. In the 
e —> 0 + limit, the joint node-layer centrality trajectory is given by aiM^ chain \ which 
we show using a dashed line. Interestingly, for the e values that exhibit slowly-varying 
conditional-node-layer-centrality trajectories in panel (a) (i.e., e < 10~ 3 ), we find that 
the joint node-layer centralities have a similar order of magnitude across the T time 
layers. For example, when e = 10~ 4 , the joint node-layer centrality of Georgia Tech at 
time t = 1945 is roughly one tenth that of Georgia Tech at time t = 1965. In contrast, 
when e > 10 _1 , the joint node-layer centralities for Georgia Tech are concentrated 
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Table 5.2 

Top time-averaged centralities (ai) and first-order-mover scores (mi) for our temporal gener¬ 
alization [see Eg. (|3.1|)/ of PageRank with p = 0.85 for universities in the Ph.D. exchange network. 


Top Time-Averaged Centralities Top First-Order-Mover Scores 


Rank 

University 


mi 

University 

rrn 

OLi 

1 

MIT 

0.484 

1.82 

Princeton 

15.01 

0.303 

2 

Stanford 

0.402 

7.58 

UC Berkeley 

10.76 

0.390 

3 

UC Berkeley 

0.390 

10.76 

Harvard 

7.99 

0.363 

4 

Harvard 

0.363 

7.99 

Stanford 

7.58 

0.402 

5 

Princeton 

0.303 

15.01 

U Michigan 

5.17 

0.101 

6 

Cornell 

0.142 

2.11 

Columbia 

5.03 

0.088 

7 

U Chicago 

0.122 

4.28 

U Chicago 

4.28 

0.122 

8 

Yale 

0.112 

2.36 

NYU 

3.28 

0.102 

9 

Wisconsin-Madison 

0.111 

2.09 

Carnegie Melon 

2.83 

0.088 

10 

NYU 

0.102 

3.28 

U Washington 

2.67 

0.053 


at just a few time layers near t = 1982. Note that the dominant eigenvalue of the 
centrality matrix CK 1982 ) for this time layer is larger than those for the other time 
layers. Thus, we find that the joint node-layer centralities to localize at layer t = 1982 
as e —> oo, a phenomenon that can be described using Perron-Frobenius theory 87 


This localization implies that a single time layer is dominating the joint node-layer 
centralities, which we confirm with the observation that the marginal layer centralities 
are also localized at layers near t = 1982 (not shown). 

5.1.3. MGP: PageRank Centrality. We also examine temporal centralities 
in the MGP network using centrality matrices given by PageRank matrices 


P (t) = pA (t) diag[dJ°, • • •, d^]" 1 + (1 - p)v 1 


(5.1) 


where d 


(t) 


= £*4 


is the out-degree of node j, the quantity l — p £ [0,1] is the 
so-called teleportation parameter (also called the damping coefficient), 1 is a vector 
of ones, and v is the personalized PageRank vector (which we set to be v = iV” 1 1). 
We handle dangling nodes (i.e., nodes with out-degree 0) by adding a single self¬ 
edge for each of these nodes. (This process of adding self-edges is sometimes called 
sink preferential PageRank 46 .) By adding these edges, A (t Miag[d^ , 


4] 1 is 


guaranteed to be column stochastic, so its spectral radius is 1. PageRank centrality is 
equal to the right dominant eigenvector of the PageRank matrix. In our usage, 
in Eq. (3.1) is given by P^, and the temporal extension of PageRank is the right 


dominant eigenvector of the corresponding supra-centrality matrix <C(e). 


In Table 5.2 we show the universities with top time-averaged PageRank centrality 
and first-order-mover scores in the MGP Ph.D. exchange network. We identify similar 


top universities with PageRank as we did with authority scores (see Table 5.1). For 
example, both tables identify MIT as the university with the top time-averaged cen¬ 
trality. Interestingly, the gap between oti for MIT and the second-ranked university is 
much smaller for PageRank than it is for authority scores. We hypothesize that one 
major factor that contributes to this difference is the impact of self-edges. In par¬ 
ticular, we identify 172 self-edges for MIT during the period 1946-2010, whereas the 
second-largest number of self-edges is only 59 (for both Stanford and UC Berkeley). 
This suggests that it may be interesting to examine variants of PageRank that use a 
different strategy to deal with dangling nodes. 
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5.1.4. MGP: Summary. Our case study illustrates practical considerations 
and techniques that are useful for understanding centrality in temporal (and other 
multilayer) networks. In our analysis of authority centralities, we identified two im¬ 
portant characteristics of centrality trajectories: “slow variation” and “layer local¬ 
ization”, which we observe, respectively, in the limits e —>■ 0 + and e —>■ oo. As we 
have demonstrated, a qualitative comparison of centrality trajectories to these limit¬ 
ing cases is help ful fo r quantifying regimes of large and small e (e.g., e > 10” 1 and 
e < 10 -3 for Fig. 5.2). In this example, we found for sufficiently small e that the cen¬ 
trality trajectories vary slowly with time. As e increases, we found a transition, which 
we observed in two ways: the joint centrality localizes to just a few layers, and the 
conditional centrality begins to exhibit large fluctuations from one time layer to the 
next (i.e., trajectories no longer vary slowly). We believe this weak-coupling regime to 
be inappropriate for the MGP Ph.D. exchange network, as mathematics department 
prestige should not fluctuate wildly from one year to the next 115 . Instead, it should 


change on a slower time scale. The strong-coupling regime is described by our sin¬ 
gular perturbation (e —> 0 + ) analysis. The transition between the weak-coupling and 
strong-coupling regimes can be rather complicated. For example, see our calculations 
for Georgia Tech in Fig. |5.2| Obtaining a complete description of the dependency 
on e of centrality trajectories for all universities (i.e., not just Georgia Tech) is even 
more complicated. For scenarios in which exploring various e is not computationally 
feasible, we highlight that restricting one’s attention to the limit e —> 0 + can still yield 
very informative results (e.g., see Fig. 5.1), and obviously it is much more efficient 
computationally. 


We have also observed fascinating phenomena, and such phenomena and our tech¬ 
niques for investigating them provide avenues for further study. One phenomenon 
is eigenvector localization for multilayer networks. Eigenvector localization is well- 
known and has received considerable attention for time-independent, single-layer net¬ 
works 17,22|[31,37][47 68f86p8 89p8| l04|. Localization can arise from various forms 
of structural heterogeneity—including the presence of high-degree nodes, community 
structure, clustering, core-periphery structure, and edge weighting. For the purpose 
of ranking nodes with an eigenvector-based centrality, localization can sometimes be 
problematic, because the centrality concentrates onto a (potentially very small) subset 
of the nodes, and this makes it difficult to reliably rank nodes outside of that subset. 
This has prompted the introduction and investigation of new centrality measures that, 
for example, do not exhibit localization (or at least exhibit less severe localization) 
due to the presence of nodes with large degree 68,86 1 , although localization can still 
arise due to other network structures 98 1 . We also remark that the “non-backtracking 
centrality” (also called “Hashimoto centrality”) introduced in Ref. 86 is based on an 
eigenvector, and the framework of the present paper thus allows us to generalize it 
for temporal networks. 


For our numerical experiments in Fig. |5.2[ in the limit of large e, we observe local¬ 
ization of the eigenvector -v(e) onto layer t = 1982, which is the layer whose centrality 
matrix Cfi 1982 ) has the largest spectral radius. Given that the layers decouple in the 
limit e —> oo, and M(e) thus becomes block diagonal in that limit, localization to 
layer t = 1982 is well-described by Perron-Frobenius theory. Interestingly, we find 
for moderate e (e.g., e = 0.01) that localization appears to first occur for layer 1967, 
which also corresponds to a centrality matrix with a large spectral radius. Therefore, 
localization in this context depends both on a supra-centrality matrix’s block-diagonal 
structure (e.g., similar to localization with community structure 17|) and on the pres- 
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Table 5.3 

Top time-averaged centralities (oti) and first-order-mover scores (mi) for actors during the 
Golden Age of Hollywood (GAH). We give results for our temporal generalization [see Eq. H3/ of 
authority scores \73f . 



Top Time-Averaged Centralities 

Top First-Order-Mover Scores 

Rank 

Actor 

Oii 

rrii 

Actor 

rrn 

Oii 

1 

Gable, Clark 

0.3683 

136.32 

Marx, Groucho 

163.34 

0.3627 

2 

Marx, Groucho 

0.3627 

163.34 

Gable, Clark 

136.32 

0.3683 

3 

Marx, Harpo 

0.2844 

112.28 

Marx, Harpo 

112.28 

0.2844 

4 

Garland, Judy 

0.2820 

100.28 

Garland, Judy 

100.28 

0.2820 

5 

Tracy, Spencer 

0.2681 

98.20 

Tracy, Spencer 

98.20 

0.2681 

6 

Stewart, James 

0.2371 

78.78 

Crawford, Joan 

90.58 

0.2369 

7 

Crawford, Joan 

0.2369 

90.58 

Marx, Chico 

86.39 

0.2055 

8 

Astaire, Fred 

0.2103 

73.29 

Stewart, James 

78.78 

0.2371 

9 

Marx, Chico 

0.2055 

86.39 

Astaire, Fred 

73.29 

0.2103 

10 

Cagney, James 

0.1779 

69.00 

Cagney, James 

69.00 

0.1779 


ence of nodes with high degr ee (which con tribute to large eigenvalues in the centrality 
matrices {C^} of the layers 106|l2l fl22 ). These preliminary findings identify eigen¬ 
vector localization in multilayer networks as an exciting direction for further study. 
We do not observe a similar localization phenomenon for our temporal generalization 
of PageRank, because the spectral radii of the PageRank centrality matrices that we 
couple into supra-centrality matrices are all equal to 1, so no single layer dominates. 
Given this observation, we expect that it can be beneficial for certain applications to 
weight the layers’ centrality matrices to control their relative importances. 

5.2. Top Billing in the Golden Age of Hollywood (GAH). In our second 
case study, we examine the centralities of actors who are known for their performances 
during the so-called “Golden Age of Hollywood” (GAH)—a time period spanning 
roughly 1920-1960. To focus on the most important actors of this period, we restrict 


our study to 55 movies stars (26 female and 29 male) identified by Wikipedia 126 


as being notable within Hollywood’s Golden Age. For these individuals, we use the 
Internet Movie Database (IMDb) [63] to define a weighted, directed, temporal network 
in which each node represents a movie star and each edge encodes the number of 
times that a pair of individuals costarred in a movie during a given time window. We 
consider all movies involving this set of movie stars, and we bin the data by decade 
over the time period 1909-2009. That is, for each lO-yeai^time window t. we include a 
directed edge (and add a unit weight to the edge) * —> j (with i ^ j) for each instance 
in which the billing position of actor j is equal to or higher than that of i (i.e., actor j 
appears earlier in the credits). If the relative billing position is unknown, we add unit 
weight to the edges mj (i.e., we include both directed edges). It would be interesting 
to explore the effect of different binning and network-construction strategies. We have 
made the temporal network data available as Supplementary Material and at |l8j. 

Because the temporal GAH network is directed, we again choose to use a supra- 
centrality matrix in which the authority matrices are along the diagonal blocks. In 
Table |5.3[ we list the actors with the top time-averaged authority centralities and 
the top first-order-mover authority scores. All of these actors appeared in some of 


5 To incorporate all available data in IMDb, our first layer represents the 11-year time window 
1909-1919. All subsequent layers correspond to 10-year time windows. 
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Fig. 5.3. (Top) Conditional-centrality and (bottom) joint-centrality trajectories for the actors 
in Table fO| with the highest time-averaged centralities. The columns depict centrality trajectories 
for various values of e, which increase from left to right. Note that although Clark Gable has the 
highest time-averaged centrality, there exist values of e (e.g., e = 0.001 ) for which Groucho Marx has 
higher conditional and joint centralities at all instances in time, indicating that there is insufficient 
evident in our study to conclude objectively which actor is more significant. 


their most famous roles during the 1930s and 1940s, and many of these roles led to 
prestigious awards (e.g., Oscar nominations and wins). Interestingly, the ten actors 
actors with highest time-averaged centralities also have the highest first-order-mover 
scores, and we thus do not identify any major shifts in the centrality trajectories of 
these top actors over time. During this time period, men were nearly always billed 
higher than women, so it is not surprising that only two women appear in Table |5.3[ 
Further study of the GAH network reveals additional insights about actor central¬ 
ities. Specifically (see Table 5.3), Clark Gable has the highest time-averaged centrality 
at (0.3683), but the value (0.3627) for Groucho Marx, and both of these values are 


much larger than that (0.2844) for the third-ranked actor, Harpo Marx. In Fig. 5.3 


we plot the conditional-centrality and joint-centrality trajectories for these three ac¬ 
tors. The columns show centralities for several choices of e, which increase from left to 
right. For e = 0.00001 (left column), the relative centralities—that is, the rank order 
of actors and the ratio of one actor’s centrality versus any other’s—are very similar 
to the limiting values (see Table [5~3] ) that one obtains for e —> 0. However, there exist 
values of e in which the centrality trajectory for Groucho Marx is consistently higher 
than that for Clark Gable. For example, see the second column (for which e = 0.001) 
in Fig. |5.3| This observation supports our belief that the examination of a single cen¬ 
trality measure (or, more generally, methodology for ranking) is insufficient to fully 
describe the importances of nodes in a network. In particular, ranking nodes based on 
any centrality calculation is sensitive to the choices that are made when defining the 
centrality. Investigating the robustness of centrality measures is an important open 
question in network science (see e.g., |36l). 


5.3. Citations of United States Supreme Court Decisions (SCD). In our 

final case study, we investigate the interconnectedness of Supreme Court decisions 
(SCD) in the Unites States by examining networks that encode citations of decisions 
|8|41|42|79| . Such an investigation can reveal a variety of insights about the decisions, 
including identifying which ones build on one another and illuminating the rise and 
fall of importance of decisions. One can also try to reveal insights into large-scale 
social processes during a given time period (e.g., the identification of which social 
issues are considered to be important and/or controversial). 
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Table 5.4 

Top time-averaged, centralities (at) and first-order-mover scores (mi) for Supreme Court deci¬ 
sions \6f\. We give results for our temporal generalization [see Eq. (|3.1|)/ of authority scores 17 4 . 



Top Time-Averaged Centralities 



Rank 

Decision 

Oii 

rrn 

1 

Gibbons v. Ogden (22 U.S. 1, 1824) 

0.172 

284.8 

2 

Minnesota Rate Case (230 U.S. 352, 1913) 

0.160 

214.2 

3 

McCulloch v. Maryland (17 U.S. 316, 1819) 

0.157 

173.3 

4 

Brown v. Maryland (25 U.S. 419, 1827) 

0.107 

212.7 

5 

Robbins v. Shelby County Tax. Distr. (120 U.S. 489, 1887) 

0.096 

192.0 

6 

Cooley v. Board of Wardens (53 U.S. 12, 1851) 

0.096 

191.2 

7 

Ex parte Young (209 U.S. 123, 1908) 

0.086 

53.8 

8 

Galveston & S.A. Ry. Co. v. Texas (210 U.S. 217, 1908) 

0.083 

126.6 

9 

Cantwell v. Connecticut (310 U.S. 296, 1940) 

0.082 

516.2 

10 

Welton v. State of Missouri (91 U.S. 275, 1875) 

0.082 

168.5 


Top First-Order-Mover Scores 


Rank 

Decision 

mi 

on 

1 

Cantwell v. Connecticut (310 U.S. 296, 1940) 

516.21 

0.0822 

2 

Schneider v. State (308 U.S. 147, 1939) 

439.73 

0.0707 

3 

Thornhill v. Alabama (310 U.S. 88, 1940) 

388.08 

0.0633 

4 

Lovell v. City of Griffin (303 U.S. 444, 1938) 

369.00 

0.0632 

5 

Near v. Minnesota (283 U.S. 697, 1931) 

344.06 

0.0729 

6 

Gitlow v. New York (268 U.S. 652, 1925) 

316.33 

0.0626 

7 

DeJonge v. Oregon (299 U.S. 353, 1937) 

310.63 

0.0554 

8 

Stromberg v. California (283 U.S. 359, 1931) 

306.68 

0.0563 

9 

Chaplinsky v. New Hampshire (315 U.S. 568, 1942) 

302.71 

0.0487 

10 

Whitney v. California (274 U.S. 357, 1927) 

291.03 

0.0582 


We study the data set that was made available by Fowler et al. 40 41 , and we 


note that the complete decisions are available online from the U.S. government 64 


We examine temporal citation networks for the time range 1800-2002, which we bin 
into decades to give T = 20 time layers^ We construct a directed temporal network 
in which we include a directed edge from node i to node j at time t if decision i cites 
decision j and decision i was written during the tth decade. To study centrality in 
such a temporal network, we restrict our attention to the largest weakly-connected 
component, which contains N = 25, 389 nodes. We study our temporal generalization 
of authority scores 73 , as high authority nodes should correspond to highly-cited, 


influential decisions. 

In Table [5~H we indicate the decisions that have the top time-averaged authority 
scores and those that have the top first-order-mover authority scores. We identify the 
top-ranking decision to be Gibbons v. Ogden (22 U.S. 1, 1824), which is well-known 
to be a highly influential commerce case 4ll 42]. (Note that we use standard case 


notation, so for this example citation, 22 is the volume, 1 is the page, and 1824 is 
the year.) The decisions with the top time-average centralities tend to be decisions 
from before 1900. In contrast, the nodes with the top first-order-mover scores tend 
to be decisions from the period 1920-1940. These decisions initially have very low 


The final time layer is slightly longer than a decade, because it encompasses the years 1990-2002. 
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centralities because they do not exist early in the network, but influential decisions 
from this period later achieve high authorities. For example, see Table 2 in Ref. 41 


which identified Cantwell v. Connecticut (310 U.S. 296, 1940) as the node with top 
authority for a network in which there exists an edge i —► j if and only if decision i cites 
decision j during the years t £ [1754, 2002]. This decision is the only one that makes 
both of our top-10 lists in Table [5~i| Interestingly, the cases with top time-averaged 
centralities deal mostly with interstate commerce and the limits of state governments 
in regulating it, and the cases with top first-order-mover scores deal primarily with 
First-Amendment rights. 


6. Conclusions. We developed and analyzed a generalization of eigenvector- 
based centrality measures for temporal networks, and we demonstrated the utility of 
such temporal centralities for identifying important entities in three case studies: the 
Ph.D. exchange network for mathematical sciences in the United States, costarring in 
the Golden Age of Hollywood, and citations of decisions in the United States Supreme 
Court. Consistent with the lessons from the development of multilayer generalizations 
of modularity |2j[5,92 , an essential ingredient of generalizing centrality measures to 
multilayer representations of temporal networks is to give different treatment to intra¬ 
layer edges between different physical nodes and inter-layer edges that connect the 
same node across time. We incorporated inter-layer edges by constructing a supra- 
centrality matrix, which requires a network to either have discrete time or be binned 
(i.e, aggregated) into discrete times using time windows (e.g., using summation or 
another process [120 123 ). Constructing this matrix requires making a choice of 
centrality matrix (e.g., authority, hub, adjacency, PageRank, and so on), and one 
then couples these matrices using an inter-layer coupling of strength u> (i.e., 1/e). 

We showed that the dominant eigenvector of a supra-centrality matrix charac¬ 
terizes the “joint centralities” of node-layer pairs, and we introduced the concepts 
of “marginal centralities” and “conditional centralities” to study (1) the decoupled 
importances of nodes and layers and (2) the nodes’ centrality trajectories over time. 
One needs to be cautious about interpreting joint centralities, as there is a bias to¬ 


wards the middle layers [e.g., see Figs. |3.1[ a) and 5.2 ]b)], but conditional centralities 
do not seem to have a similar bias. In our temporal centralities, we observe properties 
such as localization (which is an important issue for centrality measures more gen¬ 
erally [68 86,98 120| ) and different time scales for how node centrality changes over 
time. For example, we observed eigenvector localization, amounting to localization 
in time layers, in the limit w —> 0. Further research is important to explore eigen¬ 
vector localization in temporal and multilayer networks and to examine the results 
from generalizing different types of eigenvector-based centralities. It would also be 
worthwhile to explore whether concepts from statistics about smoothing, such as us¬ 
ing cross-validation to choose bandwidth parameters, can be used to help guide the 
selection of values of u. 


By focusing on the strong-coupling limit— including the construction of a per¬ 
turbation expansion in this singular limit—we derived simple, principled formulas 
to define time-averaged centralities and first-order-mover scores (which measure the 
magnitude of a physical node’s centrality changes over time). This makes it possible 
to easily identify not only which entities are most central in a temporal network but 
also which ones are the “top movers” in centrality over time. Our methodology works 
for any eigenvector-based centrality, which we define as centrality measures in which 
the nodes’ centralities are given by the entries of the dominant eigenvector of a ma¬ 
trix. There are numerous popular types of eigenvector-based centralities (including 
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PageRank centralities 95], hub and authority centralities 73 , and eigenvector cen¬ 
trality and new centralities of this form continue to be developed [ 86 ;. We have 


also provided Matlab software 118 that implements our methodology. It can be 


used to extend any eigenvector-based centrality to temporal networks. 

An important direction for future work is to explore different strategies for the 
construction of inter-layer edges. We have studied inter-layer edges between nearest- 
neighbor node-layer pairs; that is, a node-layer pair (i,f) is adjacent to (i,t — 1 ) and 
(i,t + l), so the edges form an undirected chain that bridges physical node i across the 
T time layers. Other strategies should be explored (see discussions in 721), including 
ones with inter-layer edges that are directed (e.g., other temporal generalizations 
of communicability centrality [38]). Directed inter-layer edges are able to represent 
causal scenarios, but such causal coupling can yield (supra-centrality) matrices that 
are not irreducible, which is a problematic situation for eigenvector-based centrality 
measures. Specifically, one would need to construct an associated network structure 
that satisfies strong connectivity (e.g., using ideas such as teleportation 46 , 77]). 


Developing eigenvector-based centrality measures that respect causality is thus one 
exciting direction, although one then needs to keep careful track of the biases caused by 
the choice of teleportation strategy. We note that we have developed our perturbation 
analyses for general inter-layer coupling, so they can be applied to such scenarios. 
Finally, it is worth exploring our approach for different applications such as change- 
point detection and for a variety of example networks such data streams and as those 
with a large number of time layer. 


Appendix: Higher-Order Terms. In Secs. |4.2| and |4.3 


we derived zeroth- 
for 


order and first-order solutions to the dominant-eigenvalue equation, Eq. (4.6 1 
the supra-centrality matrix. 


These led to principled expressions for time-averaged 

We now 
Such 


centralities [see Eq. (4.12)] and first-order-mover scores [see Eq. (4.26)]. 
derive and solve higher-order terms in our singular perturbation expansion, 
terms can be useful for approximating v(e) for fixed e > 0. Similar to the expressions 
that we derived in Sec. [4] we obtain expressions in the form of linear equations of 
dimension N. For longitudinal data, these equations are much more computationally 


efficient to solve than directly solving Eq. (4.6), which is an eigenvalue equation of 
dimension NT. 


From the eigenvalue equation given by Eq. (4.6), we now derive fcth-order ex- 

' k Ak and v(e) 


E k 
7 — 


j—o 


E k 
3= 0 


e k v k . 


We can consider 


pansions of the form A(e) 
arbitrarily large nonnegative integers k. We derived zeroth-order (k = 0) and first- 
order (k = 1) approximations in Sec. [4] and we now derive the second-order (fc = 2) 
expansion. Because we already showed that A 2 = a 7 '^ 2 ' 1 a in Eq. ( |4.23 l, all that is 
left is to derive an expression for v 2 . Starting from Eq. (4.17), we write 


v 2 = LjLxV! + ^ , 


( 6 . 1 ) 


where the sum with constants 7 ,- = (Puj) T v 2 accounts for the projection of v 2 onto 
the null space of Lq. We solve for the constants { 7 ^} by examing the third-order 
(k = 3) expansion, which leads to 


L 0 v 3 = Liv 2 - A 2 vi - A 3 v 0 . 


( 6 . 2 ) 
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In general, 


Lo^fc = Li^fc-i — ^2 ^j^k-j ■ 
i=2 


(6.3) 


We left-multiply Eq. (6.2 1 by (Pu,:) T and note that Pup is in the null space of L 0 to 
obtain 


0 = (Pup) T Liv 2 - A 2 Pi - A 3 cp 


(6.4) 


Using the solution for v 2 given by Eq. (6.1), it then follows that 

- y^7j(PiUi) T LiPmj = (Puj) T LiL ( ! ) Li'Vi - A 2 Pi - A 3 cp 
j 

Recalling that (Pup) T LiPuj = X-p — Ai Sij, Eq. (6.5) becomes 


(6.5) 


- (x^ - Ail) 7 = g - A 2 /3 - A 3 a , 


( 6 . 6 ) 


where q = [g 1? ..., g_/v] T and qi = (Ptu^) t L 1 LqL 1 v 1 . We now solve for A 3 . Using the 
fact that a is in the null space of (X^ 1 ) — Ail), left-multiplication of Eq. ( 6 . 6 ) by a 1 
gives 


A 3 = a 1 q = tv^LiLJLi'Vi , 


(6.7) 


where we have used the relation a T (3 = 0. As before, we solve for 7 using the 
Moore-Penrose pseudoinverse: 


7 = - (X« - Ail) f (q - A 2 /3 - A 3 a) + cq , 


( 6 . 8 ) 


where c = a T 7 is the (possibly nonzero) projection of 7 onto a. 

Similar to our calculation of b = at T (3 = 0 in Sec. 4.3 we solve for the constant c 
by examining the norm of the vector v(e) ss wq + ewi + e 2 - v 2 . We require that 


1 = 


+ evi + e 2 v 2 || 2 


(6.9) 


We substitute 
terms to obtain 


"Voll + e ”ll w i|| + 2e(wo, vi) + 2e (vq, v 2 ) + 0(e ). 

= 1 and (vo,vi) = a T (3 = 0 into Eq. (6.9) and equate the 0(e 2 ) 


c = (v 0 W 2 ) = -*lki|| 2 . 


( 6 . 10 ) 


This completes the second-order approximation. 

We now derive the third-order (k = 3) approximation. Having solved A 3 in 
Eq. (6.7), we seek to find an expression for w 3 . Starting from Eq. (6.2), we write 


w 3 = LjLiw 2 - AoLj^i + ^2 


( 6 . 11 ) 
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where the sum with constants = (Ptu.,') t v 3 accounts for the projection of w 3 onto 
the null space of Lj. Using the general form given by Eq. (6.3) for k = 4 yields 


IL0W4 — L1V3 — A2W2 — A3V1 — A4W0 • 
We left-multiply Eq. (6.12) by (Puii) T to obtain 

0 = (Pu i ) T Li’v 3 - A 2 7 i - A 3 Pi - A 4 ctj. 


( 6 . 12 ) 


(6.13) 


We next substitute our solution for v 3 given by Eq. (6.11) into Eq. (6.13) to obtain 
- ^(PtgfLiPm, = (Pu i ) T LiLj (Liv 2 - A 2 w 4 ) - A 2 y» - A 3 /3 ) - A 40 *, (6.14) 


which in matrix notation is 


(X« - A 4 I)£ = r - A 27 - A 3 /3 - A 4 c* , (6.15) 

where we define r = [r 4 ,..., and r, = (Pu i ) T L 1 LQ (L 4 v 2 — A 2 v 4 ). As before, we 
use the Moore-Penrose psuedoinverse to obtain 

| = (XW - Ail)^ - A 2 (XW - A 4 I)t 7 - A 3 (X^) - Ail)^ + da , (6.16) 

which uses the fact that a is in the null space of (X^ — A 4 I). The constant d = a T £ 
gives the projection of £ onto a , and we calculate it by constraining the third-order 
expansion for v(e) to have a norm of 1. This yields 

1 = ||^o + e^i + e 2 v 2 + e 3 -w 3 || 2 
= |Ko|| 2 + e 2 |Ki|| 2 + 2e(v 0 , v 4 ) + 2e 2 (v 0 , v 2 ) 

+ 2e 3 (v 0 ,v 3 ) + 2e 3 (^ 1 ,v 2 ) + 0(e 4 ). (6.17) 


Recall that normalization of the zeroth-, first-, and second-order expansions, respec¬ 
tively, yields the relations ||v 0 || 2 = 1 , (^oWi) = 0 , and (vo,¥ 2 ) = —(l/ 2 )||^i|| 2 . 
Equating the third-order terms in Eq. (6.17) then necessitates that 


d = (¥ 0 ,¥ 3 ) = - (vi,v 2 ). 


(6.18) 


This completes the third-order approximation. 

In Figs. 6.1 ’a) and |6.1[ b), we validate the above results using the MGP and GAH 
networks, which we studied in Sec.|5.1|and Sec.|5.2[ respectively. We plot the L 2 error 


for the fcth-order approximation ||-w(e) — JX =1 dWfc || 2 to the dominant eigenvector. 
We plot these approximations for several choices of k and several values of e. When 
e is sufficiently small (i.e., e£ 3x 10 -4 for the MGP and e ;$ 10 -3 for the GAH), 
we observe (as expected) that the error decreases with increasing k. We derived our 
approximate expressions in the limit e —X 0 + , so we only expect them to be accurate 
for sufficiently small e (although asymptotic expressions do have a long history of 
often being accurate even in many situations in which there are no guarantees for 
such success 56 ). We also obtain the expected decay rates in the error as e —> 0 + . 
We find linear decay for the zeroth-order approximation, quadratic decay for the 
first-order approximation, and so on. 
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(a) Approximation accuracy MGP (b) Approximation accuracy GAH 




Fig. 6.1. We show, for several values of inter-layer coupling weight e > 0, the accuracy of 
kth-order approximate solutions to Eq. ( |4.6[ ) for k E {0,1,2,3}. We show results for (a) the MGP 
network studied in Sec. E3 and (b) the GAH network studied in Sec. \ 5. 2\ In both panels, we 
measure the error by calculating the L 2 norm ||v(e) — . =1 ^ II 2 of the difference between the 

approximate and actual dominant eigenvector. As expected, we find for sufficiently small e that the 
error decreases with increasing approximation order. The decay rate of the L 2 error as e —> also 
follows the expected scaling. 
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